Data from 1960 to 2010 raining days, there are 4168 observations 13 days of which are rain exceeds 1000 (10cm)

Change of exceeds 1000 is 0.003119002 0.31% change

install.packages("moments")                 # Install moments package
trying URL 'https://cran.ma.imperial.ac.uk/bin/macosx/contrib/4.0/moments_0.14.tgz'
Content type 'application/x-gzip' length 54265 bytes (52 KB)
==================================================
downloaded 52 KB

The downloaded binary packages are in
    /var/folders/4x/m08khw3x0497mht_rr4wwxmh0000gn/T//Rtmpmqs2am/downloaded_packages

Data

read.csv("./raw/prcp.csv") %>%
  mutate(date = as.Date(date)) %>%
  mutate(year = format(date, "%Y")) %>%
  filter(prcp != 0) %>%
  dplyr::select(date, year, prcp) -> data

data %>%
  filter(year <= 2010) %>%
  pull(prcp) %>%
  hist(breaks = 1000, plot = FALSE) -> hist

Modling Distribution

Use data from 1960 - 2010 for modeling precipitation

hist$counts %>% sum() -> sum
x <- hist$mids 
y <- hist$counts/sum 

data.frame(x, y) -> dens 
model <- drm(data = dens, y ~ log(x), fct = EXD.2())
summary(model)

Model fitted: Exponential decay with lower limit at 0 (2 parms)

Parameter estimates:

                Estimate Std. Error t-value   p-value    
d:(Intercept) 0.15567499 0.00099172  156.97 < 2.2e-16 ***
e:(Intercept) 1.27424832 0.00578594  220.23 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error:

 0.001088457 (945 degrees of freedom)
d = 0.15567499
e = 1.27424832

k <- 1/e/d
plot(x = dens$x, y = fitted(model), type = "l", col = "blue", xlim = c(0, 300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))
mod.exd <- drm(y ~ x, fct = EXD.2(), data = dens)
plot(x = dens$x, y = fitted(mod.exd), type = "l", col = "blue", xlim = c(0,300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))
summary(mod.exd)
1 - pexp(log(1000)*k, d)

1/d

data %>%
  filter(year <= 2010) %>%
  pull(prcp) %>%
  mean() %>%
  log()*k

0.31% change versus 0.442%

Labs:

A Single Years

yr = 1960
data %>%
  filter((year >= yr) & (year <= yr + 9)) %>%
  pull(prcp) %>%
  hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids 
y <- hist$counts/sum 
data.frame(x, y) -> dens 
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)
yr = 1970
data %>%
  filter((year >= yr) & (year <= yr + 9)) %>%
  pull(prcp) %>%
  hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids 
y <- hist$counts/sum 
data.frame(x, y) -> dens 
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)
yr = 1980
data %>%
  filter((year >= yr) & (year <= yr + 9)) %>%
  pull(prcp) %>%
  hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids 
y <- hist$density 
data.frame(x, y) -> dens 
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)

Risks Through out Years

mtr <- data.frame(year = YEAR)
list <- c()

for (i in YEAR){
    data %>%
      filter((year >= i) & (year <= i + 9)) %>%
      pull(prcp) %>%
      hist(breaks = 400, plot = FALSE) -> hist
    hist$counts %>% sum() -> sum
    x <- hist$mids 
    y <- hist$counts/sum 
    data.frame(x, y) -> dens 
    mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
    d <- mod$parmMat[1]
    e <- mod$parmMat[2]
    k <- 1/e/d
    p <- 1 - pexp(log(1000)*k, d)
    list <- c(list, p)
}
plot(1:length(list), list)

Use smaller breaks

Find the right break numbers

YEAR <- c(1960, 1980, 2000)
BREAKS <- c(19, 20, 21, 22)
lenlist <- c()
kulist <- c()
for(yr in YEAR){
data %>%
      filter((year >= yr) & (year <= yr + 19)) %>%
      dplyr::select(year, prcp) %>% 
      pull(prcp) %>%
      length() -> n
  lenlist <- c(lenlist, n)
data %>%
      filter((year >= yr) & (year <= yr + 19)) %>%
      dplyr::select(year, prcp) %>% 
      pull(prcp) %>%
      kurtosis() -> ku
  kulist <- c(kulist, ku)
}
Square_root = c()
Sturges = c()
Rice = c()
Donna = c()
for(n in lenlist){
  sq = ceiling(n^0.5)
  sg = ceiling(log(n, 2)) + 1
  rc = ceiling(2*n^(1/3))
  igh = (6*(n - 2)/(n +1)/(n +3))^0.5
  dn = log(abs(1 + kulist[which(lenlist %in% n)]/igh) , 2) + log(n, 2) + 1
Square_root = c(Square_root, sq)
Sturges = c(Sturges, sg)
Rice = c(Rice, rc)
Donna = c(Donna, dn)
}
data.frame( Year = YEAR,
          Length = lenlist,
           Square_root = Square_root, 
           Sturges = Sturges,
           Rice = Rice, 
           Donna = Donna
           ) -> meth
meth

Work Flow

In a 10 years interval compute risks of precipitation higher than 10cm (prcp > 1000) by using different breaks of histogram from usgin 250 to 1000

Choosing bin have many choices actually

YEAR <- c(1960, 1980, 2000)
BREAKS <- c(19, 20, 21, 22, 30, 40, 50, 60, 70, 80, 90)
for(n in BREAKS)
  data %>%
      filter((year >= 1960) & (year <= 1969)) %>%
      pull(prcp) %>%
      hist(breaks = n)
kurtosis(data$prcp)
list <- c()
listB <- c()
Year <- c()
dlist <- c()
elist <- c()
klist <- c()
for(n in BREAKS){
  for(i in YEAR){
    data %>%
      filter((year >= i) & (year <= i + 19)) %>%
      pull(prcp) %>%
      hist(breaks = n, plot = FALSE) -> hist
    hist$counts %>% sum() -> sum
    x <- hist$mids 
    y <- hist$density 
    data.frame(x, y) -> dens 
    mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
    d <- mod$parmMat[1]
    e <- mod$parmMat[2]
    k <- 1/e/d
    p <- 1 - pexp(log(1500)*k, d)
    list <- c(list, p)
    listB <- c(listB, n)
    Year <- c(Year, i)
    dlist <- c(dlist, d)
    elist <- c(elist, e)
    klist <- c(klist, k)
  }}
mtr <- data.frame(year = Year,
                    Breaks = listB, 
                      Probability = list, 
                  parameter_d = dlist,
                  parameter_e = elist,
                  parameter_k = klist)

Visulisation

library(wesanderson)
mtr %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  ggplot(aes(x = year, y = Probability, color = Breaks)) + 
    geom_point() + 
    geom_line() + 
    scale_color_brewer(palette = "Blues") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        ) + 
    ylab("Risks of Rain Exceeds 10cm") + 
    ggtitle("Risks of Rian Exceeding 15 cm increase evey 2 decades") +
    xlab("Two Decades Time From _ ") -> g
#library(plotly)
ggplotly(g)
n too large, allowed maximum for palette Blues is 9
Returning the palette you asked for with that many colors
mtr %>%
  mutate(year = as.factor(year)) %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  ggplot(aes(y = parameter_d, x = parameter_e, color = year)) + 
    geom_jitter(alpha = 0.8) + 
    scale_color_brewer(palette = "Reds") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        ) + 
    xlab("parameter e: skewness of the curve") + 
    ylab("parameter d: height of the curve") 

NA

Parameter d influence the height of pdf, parameter e influence the lenghth of pdf… this suggest in general, precipitation distribution cure are getting flatter and

plot_ly(x= mtr$parameter_e, y= mtr$parameter_d, z=mtr$Breaks, type="scatter3d", mode="markers", color=mtr$year)
mtr %>%
  dplyr::select(year, Breaks, parameter_d, parameter_e) %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  melt(id = c("year", "Breaks"), variable.name = "parameter") %>%
  ggplot(aes(x = year, y = value)) + 
    geom_line(aes(x = year, y = value, linetype = parameter, color = Breaks)) +
    scale_color_brewer(palette = "Purples") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        )

mtr %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  ggplot(aes(x = year, y = parameter_k, color = Breaks)) + 
    geom_point() + 
    geom_line() + 
    scale_color_brewer(palette = "Blues") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        )

calculate different

list <- c()
listB <- c()
Year <- c()
dlist <- c()
elist <- c()
for(n in BREAKS){
  for(i in YEAR){
    data %>%
      filter((year >= i) & (year <= i + 9)) %>%
      pull(prcp) %>%
      hist(breaks = n, plot = FALSE) -> hist
    hist$counts %>% sum() -> sum
    x <- hist$mids 
    y <- hist$counts/sum 
    data.frame(x, y) -> dens 
    mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
    d <- mod$parmMat[1]
    e <- mod$parmMat[2]
    k <- 1/e/d
    p <- 1 - pexp(log(600)*k, d)
    list <- c(list, p)
    listB <- c(listB, n)
    Year <- c(Year, i)
    dlist <- c(dlist, d)
    elist <- c(elist, e)
  }}
mtr <- data.frame(year = Year,
                    Breaks = listB, 
                      Probability = list, 
                  parameter_d = dlist,
                  parameter_e = elist)
mtr %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  ggplot(aes(x = year, y = Probability, color = Breaks)) + 
    geom_point() + 
    geom_line() + 
    scale_color_brewer(palette = "Blues") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        )

g <- df %>%
  ggplot(aes(x = date, y = prcp)) + 
  geom_line(color = "blue") + 
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        ) + 
  geom_line(aes(x = date, y = rep(1500, length(date))), color = "lightblue") + 
  ylab("Precipitation (in tenth of mm")
ggplotly(g)
data %>%
  filter(prcp >= 1000)
mtr %>%
  filter(Breaks == 21) %>%
  left_join(meth, by = c("year" = "Year")) %>%
  select(year, Probability, Length) %>%
  mutate(days = Probability*Length)
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQoKRGF0YSBmcm9tIDE5NjAgdG8gMjAxMCByYWluaW5nIGRheXMsIHRoZXJlIGFyZSA0MTY4IG9ic2VydmF0aW9ucyAxMyBkYXlzIG9mIHdoaWNoIGFyZSByYWluIGV4Y2VlZHMgMTAwMCAoMTBjbSkKCkNoYW5nZSBvZiBleGNlZWRzIDEwMDAgaXMgMC4wMDMxMTkwMDIgMC4zMSUgY2hhbmdlCgpgYGB7cn0KCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHJlYWRyKQpsaWJyYXJ5KGRyYykKbGlicmFyeShyZXNoYXBlMikKCmluc3RhbGwucGFja2FnZXMoIm1vbWVudHMiKSAgICAgICAgICAgICAgICAKbGlicmFyeSgibW9tZW50cyIpICMgZm9yIHNrZXduZXNzIG9mIGRhdGEKYGBgCgojIERhdGEKCmBgYHtyfQpyZWFkLmNzdigiLi9yYXcvcHJjcC5jc3YiKSAlPiUKICBtdXRhdGUoZGF0ZSA9IGFzLkRhdGUoZGF0ZSkpICU+JQogIG11dGF0ZSh5ZWFyID0gZm9ybWF0KGRhdGUsICIlWSIpKSAlPiUKICBmaWx0ZXIocHJjcCAhPSAwKSAlPiUKICBkcGx5cjo6c2VsZWN0KGRhdGUsIHllYXIsIHByY3ApIC0+IGRhdGEKCmRhdGEgJT4lCiAgZmlsdGVyKHllYXIgPD0gMjAyMCkgJT4lCiAgcHVsbChwcmNwKSAlPiUKICBoaXN0KGJyZWFrcyA9IDEwMDAsIHBsb3QgPSBGQUxTRSkgLT4gaGlzdApgYGAKCiMgTW9kbGluZyBEaXN0cmlidXRpb24KClVzZSBkYXRhIGZyb20gMTk2MCAtIDIwMTAgZm9yIG1vZGVsaW5nIHByZWNpcGl0YXRpb24KCmBgYHtyfQpoaXN0JGNvdW50cyAlPiUgc3VtKCkgLT4gc3VtCnggPC0gaGlzdCRtaWRzIAp5IDwtIGhpc3QkY291bnRzL3N1bSAKCmRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKYGBgCgpgYGB7cn0KbW9kZWwgPC0gZHJtKGRhdGEgPSBkZW5zLCB5IH4gbG9nKHgpLCBmY3QgPSBFWEQuMigpKQpzdW1tYXJ5KG1vZGVsKQpkID0gMC4xNTU2NzQ5OQplID0gMS4yNzQyNDgzMgoKayA8LSAxL2UvZApgYGAKCmBgYHtyfQpwbG90KHggPSBkZW5zJHgsIHkgPSBmaXR0ZWQobW9kZWwpLCB0eXBlID0gImwiLCBjb2wgPSAiYmx1ZSIsIHhsaW0gPSBjKDAsIDMwMCkpCmxpbmVzKHggPSBkZW5zJHgsIHkgPSBkZW5zJHksIHR5cGUgPSAicCIsIGNvbCA9IGFscGhhKCJibGFjayIsIDAuNSkpCmBgYAoKYGBge3J9Cm1vZC5leGQgPC0gZHJtKHkgfiB4LCBmY3QgPSBFWEQuMigpLCBkYXRhID0gZGVucykKYGBgCgpgYGB7cn0KcGxvdCh4ID0gZGVucyR4LCB5ID0gZml0dGVkKG1vZC5leGQpLCB0eXBlID0gImwiLCBjb2wgPSAiYmx1ZSIsIHhsaW0gPSBjKDAsMzAwKSkKbGluZXMoeCA9IGRlbnMkeCwgeSA9IGRlbnMkeSwgdHlwZSA9ICJwIiwgY29sID0gYWxwaGEoImJsYWNrIiwgMC41KSkKYGBgCgpgYGB7cn0Kc3VtbWFyeShtb2QuZXhkKQoKYGBgCgpgYGB7cn0KMSAtIHBleHAobG9nKDEwMDApKmssIGQpCgoxL2QKCmRhdGEgJT4lCiAgZmlsdGVyKHllYXIgPD0gMjAxMCkgJT4lCiAgcHVsbChwcmNwKSAlPiUKICBtZWFuKCkgJT4lCiAgbG9nKCkqawpgYGAKCjAuMzElIGNoYW5nZSB2ZXJzdXMgMC40NDIlCgojIExhYnM6CgojIyBBIFNpbmdsZSBZZWFycwoKYGBge3J9CnlyID0gMTk2MApkYXRhICU+JQogIGZpbHRlcigoeWVhciA+PSB5cikgJiAoeWVhciA8PSB5ciArIDkpKSAlPiUKICBwdWxsKHByY3ApICU+JQogIGhpc3QoYnJlYWtzID0gMTAwMCwgcGxvdCA9IEZBTFNFKSAtPiBoaXN0Cmhpc3QkY291bnRzICU+JSBzdW0oKSAtPiBzdW0KeCA8LSBoaXN0JG1pZHMgCnkgPC0gaGlzdCRjb3VudHMvc3VtIApkYXRhLmZyYW1lKHgsIHkpIC0+IGRlbnMgCm1vZCA8LSBkcm0oZGF0YSA9IGRlbnMsIHkgfiBsb2coeCksIGZjdCA9IEVYRC4yKCkpCmQgPC0gbW9kJHBhcm1NYXRbMV0KZSA8LSBtb2QkcGFybU1hdFsyXQprIDwtIDEvZS9kCjEgLSBwZXhwKGxvZygxMDAwKSprLCBkKQpgYGAKCmBgYHtyfQp5ciA9IDE5NzAKZGF0YSAlPiUKICBmaWx0ZXIoKHllYXIgPj0geXIpICYgKHllYXIgPD0geXIgKyA5KSkgJT4lCiAgcHVsbChwcmNwKSAlPiUKICBoaXN0KGJyZWFrcyA9IDEwMDAsIHBsb3QgPSBGQUxTRSkgLT4gaGlzdApoaXN0JGNvdW50cyAlPiUgc3VtKCkgLT4gc3VtCnggPC0gaGlzdCRtaWRzIAp5IDwtIGhpc3QkY291bnRzL3N1bSAKZGF0YS5mcmFtZSh4LCB5KSAtPiBkZW5zIAptb2QgPC0gZHJtKGRhdGEgPSBkZW5zLCB5IH4gbG9nKHgpLCBmY3QgPSBFWEQuMigpKQpkIDwtIG1vZCRwYXJtTWF0WzFdCmUgPC0gbW9kJHBhcm1NYXRbMl0KayA8LSAxL2UvZAoxIC0gcGV4cChsb2coMTAwMCkqaywgZCkKYGBgCgpgYGB7cn0KeXIgPSAxOTgwCmRhdGEgJT4lCiAgZmlsdGVyKCh5ZWFyID49IHlyKSAmICh5ZWFyIDw9IHlyICsgOSkpICU+JQogIHB1bGwocHJjcCkgJT4lCiAgaGlzdChicmVha3MgPSAxMDAwLCBwbG90ID0gRkFMU0UpIC0+IGhpc3QKaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQp4IDwtIGhpc3QkbWlkcyAKeSA8LSBoaXN0JGRlbnNpdHkgCmRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKbW9kIDwtIGRybShkYXRhID0gZGVucywgeSB+IGxvZyh4KSwgZmN0ID0gRVhELjIoKSkKZCA8LSBtb2QkcGFybU1hdFsxXQplIDwtIG1vZCRwYXJtTWF0WzJdCmsgPC0gMS9lL2QKMSAtIHBleHAobG9nKDEwMDApKmssIGQpCmBgYAoKIyMgUmlzayBUcmVuZHMKCiMjIFJpc2tzIFRocm91Z2ggb3V0IFllYXJzCgpgYGB7cn0KbXRyIDwtIGRhdGEuZnJhbWUoeWVhciA9IFlFQVIpCmxpc3QgPC0gYygpCgpmb3IgKGkgaW4gWUVBUil7CiAgICBkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0gaSkgJiAoeWVhciA8PSBpICsgOSkpICU+JQogICAgICBwdWxsKHByY3ApICU+JQogICAgICBoaXN0KGJyZWFrcyA9IDQwMCwgcGxvdCA9IEZBTFNFKSAtPiBoaXN0CiAgICBoaXN0JGNvdW50cyAlPiUgc3VtKCkgLT4gc3VtCiAgICB4IDwtIGhpc3QkbWlkcyAKICAgIHkgPC0gaGlzdCRjb3VudHMvc3VtIAogICAgZGF0YS5mcmFtZSh4LCB5KSAtPiBkZW5zIAogICAgbW9kIDwtIGRybShkYXRhID0gZGVucywgeSB+IGxvZyh4KSwgZmN0ID0gRVhELjIoKSkKICAgIGQgPC0gbW9kJHBhcm1NYXRbMV0KICAgIGUgPC0gbW9kJHBhcm1NYXRbMl0KICAgIGsgPC0gMS9lL2QKICAgIHAgPC0gMSAtIHBleHAobG9nKDEwMDApKmssIGQpCiAgICBsaXN0IDwtIGMobGlzdCwgcCkKfQpwbG90KDE6bGVuZ3RoKGxpc3QpLCBsaXN0KQpgYGAKCmBgYHtyfQoKYGBgCgojIyBVc2Ugc21hbGxlciBicmVha3MKCkZpbmQgdGhlIHJpZ2h0IGJyZWFrIG51bWJlcnMKCmBgYHtyfQpZRUFSIDwtIGMoMTk2MCwgMTk4MCwgMjAwMCkKQlJFQUtTIDwtIGMoMTksIDIwLCAyMSwgMjIpCmBgYAoKYGBge3J9Cmxlbmxpc3QgPC0gYygpCmt1bGlzdCA8LSBjKCkKZm9yKHlyIGluIFlFQVIpewpkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0geXIpICYgKHllYXIgPD0geXIgKyAxOSkpICU+JQogICAgICBkcGx5cjo6c2VsZWN0KHllYXIsIHByY3ApICU+JSAKICAgICAgcHVsbChwcmNwKSAlPiUKICAgICAgbGVuZ3RoKCkgLT4gbgogIGxlbmxpc3QgPC0gYyhsZW5saXN0LCBuKQpkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0geXIpICYgKHllYXIgPD0geXIgKyAxOSkpICU+JQogICAgICBkcGx5cjo6c2VsZWN0KHllYXIsIHByY3ApICU+JSAKICAgICAgcHVsbChwcmNwKSAlPiUKICAgICAga3VydG9zaXMoKSAtPiBrdQogIGt1bGlzdCA8LSBjKGt1bGlzdCwga3UpCn0KU3F1YXJlX3Jvb3QgPSBjKCkKU3R1cmdlcyA9IGMoKQpSaWNlID0gYygpCkRvbm5hID0gYygpCmZvcihuIGluIGxlbmxpc3QpewogIHNxID0gY2VpbGluZyhuXjAuNSkKICBzZyA9IGNlaWxpbmcobG9nKG4sIDIpKSArIDEKICByYyA9IGNlaWxpbmcoMipuXigxLzMpKQogIGlnaCA9ICg2KihuIC0gMikvKG4gKzEpLyhuICszKSleMC41CiAgZG4gPSBsb2coYWJzKDEgKyBrdWxpc3Rbd2hpY2gobGVubGlzdCAlaW4lIG4pXS9pZ2gpICwgMikgKyBsb2cobiwgMikgKyAxClNxdWFyZV9yb290ID0gYyhTcXVhcmVfcm9vdCwgc3EpClN0dXJnZXMgPSBjKFN0dXJnZXMsIHNnKQpSaWNlID0gYyhSaWNlLCByYykKRG9ubmEgPSBjKERvbm5hLCBkbikKfQpkYXRhLmZyYW1lKCBZZWFyID0gWUVBUiwKICAgICAgICAgIExlbmd0aCA9IGxlbmxpc3QsCiAgICAgICAgICAgU3F1YXJlX3Jvb3QgPSBTcXVhcmVfcm9vdCwgCiAgICAgICAgICAgU3R1cmdlcyA9IFN0dXJnZXMsCiAgICAgICAgICAgUmljZSA9IFJpY2UsIAogICAgICAgICAgIERvbm5hID0gRG9ubmEKICAgICAgICAgICApIC0+IG1ldGgKYGBgCgpgYGB7cn0KbWV0aApgYGAKCiMgV29yayBGbG93CgpJbiBhIDEwIHllYXJzIGludGVydmFsIGNvbXB1dGUgcmlza3Mgb2YgcHJlY2lwaXRhdGlvbiBoaWdoZXIgdGhhbiAxMGNtIChwcmNwIFw+IDEwMDApIGJ5IHVzaW5nIGRpZmZlcmVudCBicmVha3Mgb2YgaGlzdG9ncmFtIGZyb20gdXNnaW4gMjUwIHRvIDEwMDAKCkNob29zaW5nIGJpbiBoYXZlIG1hbnkgY2hvaWNlcyBhY3R1YWxseQoKYGBge3J9CllFQVIgPC0gYygxOTYwLCAxOTgwLCAyMDAwKQpCUkVBS1MgPC0gYygxOSwgMjAsIDIxLCAyMiwgMzAsIDQwLCA1MCwgNjAsIDcwLCA4MCwgOTApCmBgYAoKYGBge3J9CmZvcihuIGluIEJSRUFLUykKICBkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0gMTk2MCkgJiAoeWVhciA8PSAxOTY5KSkgJT4lCiAgICAgIHB1bGwocHJjcCkgJT4lCiAgICAgIGhpc3QoYnJlYWtzID0gbikKYGBgCgpgYGB7cn0Ka3VydG9zaXMoZGF0YSRwcmNwKQpgYGAKCmBgYHtyfQpsaXN0IDwtIGMoKQpsaXN0QiA8LSBjKCkKWWVhciA8LSBjKCkKZGxpc3QgPC0gYygpCmVsaXN0IDwtIGMoKQprbGlzdCA8LSBjKCkKZm9yKG4gaW4gQlJFQUtTKXsKICBmb3IoaSBpbiBZRUFSKXsKICAgIGRhdGEgJT4lCiAgICAgIGZpbHRlcigoeWVhciA+PSBpKSAmICh5ZWFyIDw9IGkgKyAxOSkpICU+JQogICAgICBwdWxsKHByY3ApICU+JQogICAgICBoaXN0KGJyZWFrcyA9IG4sIHBsb3QgPSBGQUxTRSkgLT4gaGlzdAogICAgaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQogICAgeCA8LSBoaXN0JG1pZHMgCiAgICB5IDwtIGhpc3QkZGVuc2l0eSAKICAgIGRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKICAgIG1vZCA8LSBkcm0oZGF0YSA9IGRlbnMsIHkgfiBsb2coeCksIGZjdCA9IEVYRC4yKCkpCiAgICBkIDwtIG1vZCRwYXJtTWF0WzFdCiAgICBlIDwtIG1vZCRwYXJtTWF0WzJdCiAgICBrIDwtIDEvZS9kCiAgICBwIDwtIDEgLSBwZXhwKGxvZygxNTAwKSprLCBkKQogICAgbGlzdCA8LSBjKGxpc3QsIHApCiAgICBsaXN0QiA8LSBjKGxpc3RCLCBuKQogICAgWWVhciA8LSBjKFllYXIsIGkpCiAgICBkbGlzdCA8LSBjKGRsaXN0LCBkKQogICAgZWxpc3QgPC0gYyhlbGlzdCwgZSkKICAgIGtsaXN0IDwtIGMoa2xpc3QsIGspCiAgfX0KbXRyIDwtIGRhdGEuZnJhbWUoeWVhciA9IFllYXIsCiAgICAgICAgICAgICAgICAgICAgQnJlYWtzID0gbGlzdEIsIAogICAgICAgICAgICAgICAgICAgICAgUHJvYmFiaWxpdHkgPSBsaXN0LCAKICAgICAgICAgICAgICAgICAgcGFyYW1ldGVyX2QgPSBkbGlzdCwKICAgICAgICAgICAgICAgICAgcGFyYW1ldGVyX2UgPSBlbGlzdCwKICAgICAgICAgICAgICAgICAgcGFyYW1ldGVyX2sgPSBrbGlzdCkKYGBgCgojIFZpc3VsaXNhdGlvbgoKYGBge3J9CmxpYnJhcnkod2VzYW5kZXJzb24pCm10ciAlPiUKICBtdXRhdGUoQnJlYWtzID0gYXMuZmFjdG9yKEJyZWFrcykpICU+JQogIGdncGxvdChhZXMoeCA9IHllYXIsIHkgPSBQcm9iYWJpbGl0eSwgY29sb3IgPSBCcmVha3MpKSArIAogICAgZ2VvbV9wb2ludCgpICsgCiAgICBnZW9tX2xpbmUoKSArIAogICAgc2NhbGVfY29sb3JfYnJld2VyKHBhbGV0dGUgPSAiQmx1ZXMiKSArCiAgdGhlbWUocGFuZWwuYm9yZGVyID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5tYWpvciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWlub3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KGZpbGwgPSAiI2ZhZmFmYSIpLAogICAgICAgICAgcGxvdC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KGZpbGwgPSAid2hpdGUiKQogICAgICAgICkgKyAKICAgIHlsYWIoIlJpc2tzIG9mIFJhaW4gRXhjZWVkcyAxMGNtIikgKyAKICAgIGdndGl0bGUoIlJpc2tzIG9mIFJpYW4gRXhjZWVkaW5nIDE1IGNtIGluY3JlYXNlIGV2ZXkgMiBkZWNhZGVzIikgKwogICAgeGxhYigiVHdvIERlY2FkZXMgVGltZSBGcm9tIF8gIikgLT4gZwojbGlicmFyeShwbG90bHkpCmdncGxvdGx5KGcpCmBgYAoKYGBge3J9Cm10ciAlPiUKICBtdXRhdGUoeWVhciA9IGFzLmZhY3Rvcih5ZWFyKSkgJT4lCiAgbXV0YXRlKEJyZWFrcyA9IGFzLmZhY3RvcihCcmVha3MpKSAlPiUKICBnZ3Bsb3QoYWVzKHkgPSBwYXJhbWV0ZXJfZCwgeCA9IHBhcmFtZXRlcl9lLCBjb2xvciA9IHllYXIpKSArIAogICAgZ2VvbV9qaXR0ZXIoYWxwaGEgPSAwLjgpICsgCiAgICBzY2FsZV9jb2xvcl9icmV3ZXIocGFsZXR0ZSA9ICJSZWRzIikgKwogIHRoZW1lKHBhbmVsLmJvcmRlciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1pbm9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIiNmYWZhZmEiKSwKICAgICAgICAgIHBsb3QuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIndoaXRlIikKICAgICAgICApICsgCiAgICB4bGFiKCJwYXJhbWV0ZXIgZTogc2tld25lc3Mgb2YgdGhlIGN1cnZlIikgKyAKICAgIHlsYWIoInBhcmFtZXRlciBkOiBoZWlnaHQgb2YgdGhlIGN1cnZlIikgCiAgICAKYGBgCgpQYXJhbWV0ZXIgZCBpbmZsdWVuY2UgdGhlIGhlaWdodCBvZiBwZGYsIHBhcmFtZXRlciBlIGluZmx1ZW5jZSB0aGUgbGVuZ2h0aCBvZiBwZGYuLi4gdGhpcyBzdWdnZXN0IGluIGdlbmVyYWwsIHByZWNpcGl0YXRpb24gZGlzdHJpYnV0aW9uIGN1cmUgYXJlIGdldHRpbmcgZmxhdHRlciBhbmQKCmBgYHtyfQpwbG90X2x5KHg9IG10ciRwYXJhbWV0ZXJfZSwgeT0gbXRyJHBhcmFtZXRlcl9kLCB6PW10ciRCcmVha3MsIHR5cGU9InNjYXR0ZXIzZCIsIG1vZGU9Im1hcmtlcnMiLCBjb2xvcj1tdHIkeWVhcikKYGBgCgpgYGB7cn0KbXRyICU+JQogIGRwbHlyOjpzZWxlY3QoeWVhciwgQnJlYWtzLCBwYXJhbWV0ZXJfZCwgcGFyYW1ldGVyX2UpICU+JQogIG11dGF0ZShCcmVha3MgPSBhcy5mYWN0b3IoQnJlYWtzKSkgJT4lCiAgbWVsdChpZCA9IGMoInllYXIiLCAiQnJlYWtzIiksIHZhcmlhYmxlLm5hbWUgPSAicGFyYW1ldGVyIikgJT4lCiAgZ2dwbG90KGFlcyh4ID0geWVhciwgeSA9IHZhbHVlKSkgKyAKICAgIGdlb21fbGluZShhZXMoeCA9IHllYXIsIHkgPSB2YWx1ZSwgbGluZXR5cGUgPSBwYXJhbWV0ZXIsIGNvbG9yID0gQnJlYWtzKSkgKwogICAgc2NhbGVfY29sb3JfYnJld2VyKHBhbGV0dGUgPSAiUHVycGxlcyIpICsKICB0aGVtZShwYW5lbC5ib3JkZXIgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1ham9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICIjZmFmYWZhIiksCiAgICAgICAgICBwbG90LmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICJ3aGl0ZSIpCiAgICAgICAgKQpgYGAKCmBgYHtyfQptdHIgJT4lCiAgbXV0YXRlKEJyZWFrcyA9IGFzLmZhY3RvcihCcmVha3MpKSAlPiUKICBnZ3Bsb3QoYWVzKHggPSB5ZWFyLCB5ID0gcGFyYW1ldGVyX2ssIGNvbG9yID0gQnJlYWtzKSkgKyAKICAgIGdlb21fcG9pbnQoKSArIAogICAgZ2VvbV9saW5lKCkgKyAKICAgIHNjYWxlX2NvbG9yX2JyZXdlcihwYWxldHRlID0gIkJsdWVzIikgKwogIHRoZW1lKHBhbmVsLmJvcmRlciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1pbm9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIiNmYWZhZmEiKSwKICAgICAgICAgIHBsb3QuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIndoaXRlIikKICAgICAgICApCmBgYAoKY2FsY3VsYXRlIGRpZmZlcmVudAoKYGBge3J9Cmxpc3QgPC0gYygpCmxpc3RCIDwtIGMoKQpZZWFyIDwtIGMoKQpkbGlzdCA8LSBjKCkKZWxpc3QgPC0gYygpCmZvcihuIGluIEJSRUFLUyl7CiAgZm9yKGkgaW4gWUVBUil7CiAgICBkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0gaSkgJiAoeWVhciA8PSBpICsgOSkpICU+JQogICAgICBwdWxsKHByY3ApICU+JQogICAgICBoaXN0KGJyZWFrcyA9IG4sIHBsb3QgPSBGQUxTRSkgLT4gaGlzdAogICAgaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQogICAgeCA8LSBoaXN0JG1pZHMgCiAgICB5IDwtIGhpc3QkY291bnRzL3N1bSAKICAgIGRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKICAgIG1vZCA8LSBkcm0oZGF0YSA9IGRlbnMsIHkgfiBsb2coeCksIGZjdCA9IEVYRC4yKCkpCiAgICBkIDwtIG1vZCRwYXJtTWF0WzFdCiAgICBlIDwtIG1vZCRwYXJtTWF0WzJdCiAgICBrIDwtIDEvZS9kCiAgICBwIDwtIDEgLSBwZXhwKGxvZyg2MDApKmssIGQpCiAgICBsaXN0IDwtIGMobGlzdCwgcCkKICAgIGxpc3RCIDwtIGMobGlzdEIsIG4pCiAgICBZZWFyIDwtIGMoWWVhciwgaSkKICAgIGRsaXN0IDwtIGMoZGxpc3QsIGQpCiAgICBlbGlzdCA8LSBjKGVsaXN0LCBlKQogIH19Cm10ciA8LSBkYXRhLmZyYW1lKHllYXIgPSBZZWFyLAogICAgICAgICAgICAgICAgICAgIEJyZWFrcyA9IGxpc3RCLCAKICAgICAgICAgICAgICAgICAgICAgIFByb2JhYmlsaXR5ID0gbGlzdCwgCiAgICAgICAgICAgICAgICAgIHBhcmFtZXRlcl9kID0gZGxpc3QsCiAgICAgICAgICAgICAgICAgIHBhcmFtZXRlcl9lID0gZWxpc3QpCgpgYGAKCmBgYHtyfQptdHIgJT4lCiAgbXV0YXRlKEJyZWFrcyA9IGFzLmZhY3RvcihCcmVha3MpKSAlPiUKICBnZ3Bsb3QoYWVzKHggPSB5ZWFyLCB5ID0gUHJvYmFiaWxpdHksIGNvbG9yID0gQnJlYWtzKSkgKyAKICAgIGdlb21fcG9pbnQoKSArIAogICAgZ2VvbV9saW5lKCkgKyAKICAgIHNjYWxlX2NvbG9yX2JyZXdlcihwYWxldHRlID0gIkJsdWVzIikgKwogIHRoZW1lKHBhbmVsLmJvcmRlciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1pbm9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIiNmYWZhZmEiKSwKICAgICAgICAgIHBsb3QuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIndoaXRlIikKICAgICAgICApCmBgYAoKYGBge3J9CmcgPC0gZGYgJT4lCiAgZ2dwbG90KGFlcyh4ID0gZGF0ZSwgeSA9IHByY3ApKSArIAogIGdlb21fbGluZShjb2xvciA9ICJibHVlIikgKyAKICB0aGVtZShwYW5lbC5ib3JkZXIgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1ham9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICIjZmFmYWZhIiksCiAgICAgICAgICBwbG90LmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICJ3aGl0ZSIpCiAgICAgICAgKSArIAogIGdlb21fbGluZShhZXMoeCA9IGRhdGUsIHkgPSByZXAoMTUwMCwgbGVuZ3RoKGRhdGUpKSksIGNvbG9yID0gImxpZ2h0Ymx1ZSIpICsgCiAgeWxhYigiUHJlY2lwaXRhdGlvbiAoaW4gdGVudGggb2YgbW0iKQpnZ3Bsb3RseShnKQpgYGAKCmBgYHtyfQpkYXRhICU+JQogIGZpbHRlcihwcmNwID49IDEwMDApCmBgYAoKYGBge3J9Cm10ciAlPiUKICBmaWx0ZXIoQnJlYWtzID09IDIxKSAlPiUKICBsZWZ0X2pvaW4obWV0aCwgYnkgPSBjKCJ5ZWFyIiA9ICJZZWFyIikpICU+JQogIHNlbGVjdCh5ZWFyLCBQcm9iYWJpbGl0eSwgTGVuZ3RoKSAlPiUKICBtdXRhdGUoZGF5cyA9IFByb2JhYmlsaXR5Kkxlbmd0aCkKYGBgCg==